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Three pseudospectral algorithms are described (Euler, leapfrog and trapez) for solving numerically 
the time dependent nonlinear Schrodinger equation in one, two or three dimensions. Numerical 
stability regions in the parameter space are determined for the cubic nonlinearity, which can be 
easily extended to other nonlinearities. For the first two algorithms, maximal timesteps for stability 
are calculated in terms of the maximal Fourier harmonics admitted by the spectral method used to 
calculate space derivatives. The formulas are directly applicable if the discrete Fourier transform 
is used, i.e. for periodic boundary conditions. These formulas were used in the relevant numerical 
programs developed in our group. 



I. INTRODUCTION 

The nonlinear Schrodinger (NLS) equation in one, two 
or three dimensions is a commonly used model in various 
branches of physics, e.g. in plasma theory 0, optics Q 
or condensed matter theory |3( • In these applications the 
cubic form of the nonlinear term is used, and our analysis 
will be pertinent to this form. However, the generaliza- 
tion to other nonlinearities will usually be straighhtfor- 
ward. To make the analysis directly applicable to various 
situations we assume the following form of the NLS equa- 
tion: 



will be linearly related to each other so as to make the 
NLS equation in machine units as simple as possible. As 
the Fourier transform of a periodic function with period 
P is the simplest for P = 2ir, we first transform the space 
intervals in J2J into the intervals [0, 2tt] by putting 



x = a x (x + L s ), 

V = Uy(y + Ly), 

z = a z (z + L z ), 
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D\u\ 2 u + Eu = Q, (1) 



where t denotes time, and x, y, z are space coordinates; 
subscripts of u denote partial derivatives. The coeffi- 
cients A (7^ 0), B, C, D 0) are assumed to be real, 
but E can be complex, E = E r + iE\. The unknown 
function will be denned in a box: 



~L X < x < L x , — Ly < y < Ly 



-L- z <z<L 2 



(2) 



with periodic boundary conditions: u(—L x , . . .) = 
u(L x , . . .), etc. These conditions will be fulfilled exactly 
for solutions periodic in x, y and z but approximately 
also for solitary solutions, which are exponentially small 
at the boundaries. Non-periodic boundary conditions re- 
quire special treatment, see the end of Sec. III. 



Furthermore we normalize the time and the unknown 
function (at, a u > 0) 



t = a t t , u — a u u . (4) 
With these transformations Eq. becomes 
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We can choose at and a u so that \A\a x /at — 1 and 
\D\/(af l a t ) = 1, leading to 



a t = \A\a 2 x , 



a,. 



V\D/A\ 



(G) 



II. TRANSFORMATION TO MACHINE UNITS 

For the machine variables, which will be used in pro- 
gramming, we will drop the bars used in to denote 
the original (physical) variables. Both types of variables 



With this choice Eq. JSJ can be written 
u t = F[u] , 

where 



(7) 



F[u] = i [sgn(A)(u xx + c y u vy + c z u zz ) 
1askor@fuw.edu.pl| +Sgn(L>) \u\ 2 u + Eu] , (8) 
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This formula defines a two-point algorithm: u(t + At) is 
given in terms of u{t) and u(t + At/2), i.e. at the centre 
of the interval [t, t + At}. Both u(t) and u(t + At/2) must 
(9) be known, and hence the actual evolution interval is At/2 
rather than At. Replacing in l|17|) u(t) — » w(t — At), the 
(10) central point will now be at t, and the integration interval 
will be 2At: 



III. NUMERICAL ALGORITHMS 

Integrating the evolution equation Q from some value 
of t to a later instant t + At, At > 0, we obtain 



u(t + At) = u(t) + 



t+At 



F[u] 



dr. 



(11) 



Assuming that At is small, simple numerical algorithms 
and their error estimates can be obtained from (|llf> by 
using the parabolic interpolation formula for F[it]| r : 

F[u]| T = F + (5{t -t)+ 7 (r - t)[r - (t + At)] , (12) 

where the quantities (3 and 7 can be expressed in terms 
of 



F = F[u] 



F 1 = F[u}\ t +At, F 1/2 = F[u]\ t+At/2 , 

(13) 



Fi — Fq 
At 
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[F - 2F 1/2 + . (14) 



(3 and 7 are nearly independent of At (close to their limits 
as At -> 0). 

Inserting (|12|l into (|llfl we can find two algorithms, 
an Euler algorithm (first order in At): 

u{t + At) = u(t) + AtF + 0[(At) 2 } , (15) 

a Trapez algorithm (second order in At): 

u(t + At) = u(t) + At i(F + Fi) + 0[(At) 3 } . (16) 

Note that JTHJl is an implicit algorithm: u(t + At) is de- 
fined in terms of u(t) and u(t + At), i.e. equation l|16|l 
must be solved for u(t + At). Usually, this is done by an 
iterative procedure, where in the lowest approximation 
Fi on the RHS of (fH))l is replaced by F . This defines 
the first approximation to u(t + At) to be used on the 
RHS of 11 6|l to define second approximation to u(t + At), 
etc. 

Another simple (explicit) algorithm can be obtained 
from (O if we put t = t + At/2 in fT2]l: 

F 1/2 = ±(F + F 1 )- 1 (At) 2 /4, 
and use h(Fo +F\) calculated from this equation in ptijl: 



u(t + At) = u(t - At) + 2AtF + 0[(At) 3 ] . (18) 

Here u(t + At) is defined in terms of u(t — At) and u(t). 
Using the known value u(t), and it (t+At) calculated from 
(|18fl we can determine u(t + 2 At), etc. (the leapfrog pro- 
cedure). The only problem is to start this procedure, 
which requires two values of u: it (to) and it (to + At), 
where to is an initial value of t. But if u(to) is pre- 
scribed, the evolution equation Q (of first order in t) 
defines u(t) for any t > to. In particular, u(to + At) is 
defined. It can be calculated up to 0[(At) 2 ] by using the 
Euler algorithm (|15|l . Using this approximation to deter- 
mine F [u] required for (|18l) . the error in F[u] will also be 
0[(At) 2 ]. After multiplication by At it will produce an 
error comparable with that in (|18l) . 0[(At) 3 ]. 

In the pseudospectral method described in [4(, the 
Discrete Fourier Transform (DFT) with respect to each 
space variable is used to calculate the derivatives in JSJ. 
Thus the interval [0, 2ir] for x will be divided into N x 
subintervals of length Ax — 2ir/N x , and similarly for y 
and z (N y subintervals of length Ay = 2tt/N v and N z 
subintervals of length Az — 2ir/N z ; the numbers N x , N y 
and N z can be either even or odd, N x — 2M X or N x — 
2M X + 1, etc.). The function u defined on the discrete 
mesh {xj,y m , z n ), Xj = jAx, y m = mAy, z n = nAz, can 
be transformed to discrete Fourier space for x, y, and z 
variables. Thus for each y m and z n we define the Discrete 
Fourier Transform in x: 



v(k x ) = -= ^ u(xj) exp(-ik x o 



(19) 



The inverse transform is given by 



u(Xj) 




v(k x ) exp(ik x x) 



(20) 



u(t + At) = u(t) + AtF 1/2 + 0[(At) 3 



(17) 



both if N x = 2M X or N x = 2M X + 1. In the first 
case, the summation index in l|20|l actually ends up with 
k x = M x — 1. However, as in that case v{k x ) exp{ik x Xj) 
is periodic as function of k x with period 2M X , Eq. I|20|) 
will be correct if only one half of the contributions at 
k x = ±M X are included in the sum over k x . Replac- 
ing x — > y, z and j — > m,n everywhere in 119fl and 120|) 
we obtain the formulas for the DFT in y (for each Xj 
and z n ) and in z (for each Xj and y m ). The essence 
of the pseudospectral approach is to calculate the partial 
derivatives at the mesh points by differentiating the inter- 
polation formula l|20|) with respect to x (or its analogues 
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with respect to y or z). Thus, for example 

7) = -7= y^(^x) 2 ^(fcx) exp(ik x Xj) , (21) 



^xx 



etc. In our numerical programs, the sums of the type l|19|) 
or (|21l) were determined by using either the Fast Fourier 
Transform (FFT) subroutine in complex domain, given 
in [j| (in the early version of the program from 1996), or 
more efficient "Multiple ID FFT subroutines for complex 
data" , taken from the NAG Fortran Library (in the later 
versions from 1999 and 2004). 

The Discrete Fourier Transform as described above is 
an efficient tool for numerical differentiation. The only 
problem is that Eq. 12U|) implies periodicity of the func- 
tion u{x) which may not be the case, neither exactly 
nor approximately. If that happens, one can follow the 
suggestion of P. J. Roach [(| to split the function u(x) 
into the periodic component given by l|2U|l and a non- 
periodic one in the form of a polynomial. Another ap- 
proach could be to replace the DFT by expansion of u(x) 
in a non-periodic orthogonal basis, e.g. that of orthogo- 
nal polynomials. This type of approach using the Cheby- 
shev polynomials has recently been discused in detail and 
tested by A. Deloff 0- 



IV. NUMERICAL STABILITY 

To examine the numerical stability of the algorithms 
derived in the previous section, we linearize F[u], Eq. JSJl, 
by replacing 



\u\ 2 u — > |uo| 2 w, wo = const . 



(22) 



As Eqs. (|15[l, (|To)) and (fT5f) (without error terms) have 
constant coefficients, their solutions can be looked for in 
the form of exponential functions of x, y, z and t: 



u = k*/ a * exp[i(k x x + k y y + k z z) 



(23) 



Numerical stability of the algorithm in question means 
that the solution 12MI) cannot grow in time, i.e., \k\ < 1. 
Inserting (|23|l into © we obtain 



F[u] = i[w + iEi\u, 

-sgn(D)\u \ 2 + E, 



w = -sgn(A)(fc 2 + Cy k 2 y + c z k 2 z ) 



(24) 
(25) 



and furthermore 

u(t + At) =uk, u(t- At) =u/k. (26) 

Eq. 123) leads to 

Mmax < M 2 X + \c y \M 2 + \c z \M 2 + \u Q \ 2 + \E r \, (27) 

where M x = max|fc x |, etc., see Eq. (|20f) . This estimate is 
an accurare expression for |w| max if all terms in l|25|l have 



the same sign. Otherwise certain terms in l)27|) should be 
discarded. Nevertheless, in practice the overestimation 
given by the RHS of 127|) in that case is not large, and if 
this formula is used in the expressions for maximal At for 
stability given in what follows, the only effect will be the 
introduction of a small safety margin. To obtain a precise 
expression for |w| max , all terms in (|25|l should be divided 
into two groups, of positive and negative. The contribu- 
tion to the RHS of i(2"Y|) of these two groups should be 
compared, and the group of terms with smaller contribu- 
tion should be discarded. 



A. Euler algorithm 

Using (|231 and J2HJ) in we obtain 

k = 1 + iAt[w + iEi] , (28) 



i.e. 



\k\ 2 = 1 - 2AtEi + (At) 2 [E 2 + w 2 } 

= 1 - At[2Ei - At{E 2 + w 2 )] . (29) 

Thus if E\ < 0, we obtain \n\ 2 > 1, i.e. the Euler algo- 
rithm is numerically unstable. 

Numerical stability is only possible for E{ > if 

< At[2Ei - At(E 2 + w 2 )} < 1 . 

As in practice |u>| ma x 3> \E{\, the stability condition for 
the Euler algorithm takes the form 



Ei > and At < 



2Ei 



E 2 + (M max ) 2 



(30) 



where |u>| max is given by (|27l) (with a possible modifica- 
tion as described above). 



B. Implicit algorithm 



Using (|33J and J23) in JTBJ) we obtain 



1 - §At(-E; 



IW) 



1 + |At(£i - iw) ' 



i.e. 



2AtE { 



(1 



\AtE,Y 



aAtwf 



(31) 



(32) 



Thus if Ei < 0, we obtain \k\ 2 > 1, i.e. the implicit 
algorithm is numerically unstable, and for E[ — this 
algorithm is marginally stable (|k| 2 — 1). And finally, 
for Ei > 0, p is positive and should not be greater than 
one (again due to expected |w| max 3> l^il), which means 
numerical stability for any value of At. 
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C. Leapfrog algorithm 
Using lO and pp ]) in l[T5 )l we obtain a quadratic in k 



2At(£i - ky)/t -1 = 0. 



(33) 



Solving Eq. I|33|) we obtain 

k = -At(Ei - iw) ± ^\ + {At) 2 {E { -iw) 2 . (34) 

For Ei ^ 0, the general expressions for the real and imag- 
inary part of k are a bit complicated, but in the limit 
At — > we easily find 



k — zbl — Ai(£; ; + iw) , i.e. 



l±2A^i. (35) 



This can always be greater than one, which means insta- 
bility. 

For E[ = 0, we obtain 



« = — ia± \/ \ — a 2 . 



a = At', 



(36) 



Hence if \a\ < 1 we obtain |k| 2 = 1, which means 
marginal stability, whereas for \a\ > 1 we get 



Nmax = |o| + \/a? - 1 > \a\ > 1 , 

which means instability. Hence the numerical stability 
condition is \a\ < 1. This condition leads to the following 
formula for maximal timestep At for stability: 



At = 



(37) 



where < c < 1, and |u>| ma x is given by (|27|l (again with 
the possible modification). 
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